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ABSTRACT 

Pulsar activity and its related radiation mechanism are usually explained by invoking 
some plasma processes occurring inside the magnetosphere, be it polar caps, outer/slot 
gaps or in the transition region between the quasi-static magnetic dipole regime and 
the wave zone, like the striped wind. Despite many detailed local investigations, the 
global electrodynamics around those neutron stars remains poorly described with only 
little quantitative studies on the largest scales, i.e. of several light-cylinder radii 7%. 
Better understanding of these compact objects requires a deep and accurate knowledge 
of their immediate electromagnetic surrounding within the magnetosphere and its link 
to the relativistic pulsar wind. This is compulsory to make any reliable predictions 
about the whole electric circuit, the energy losses, the sites of particle acceleration and 
the possibly associated emission mechanisms. 

The aim of this work is to present accurate solutions to the nearly stationary force- 
free pulsar magnetosphere and its link to the striped wind, for various spin periods and 
arbitrary inclination. To this end, the time-dependent Maxwell equations are solved in 
spherical geometry in the force-free approximation using a vector spherical harmonic- 
expansion of the electromagnetic field. An exact analytical enforcement of the diver- 
genceless of the magnetic part is obtained by a projection method. Special care has 
been given to design an algorithm able to look deeply into the magnetosphere with 
physically realistic ratios of stellar i?* to light-cylinder tl radius. However, currently 
available computational resources allows us only to set i?*/rL = 10 _1 corresponding 
to pulsars with period 2 ms. The spherical geometry permits a proper and mathe- 
matically well-posed imposition of self-consistent physical boundary conditions on the 
stellar crust. We checked our code against several analytical solutions, like the Deutsch 
vacuum rotator solution and the Michel monopole field. We also retrieve energy losses 
comparable to the magneto-dipole radiation formula and consistent with previous sim- 
ilar works. Finally, for arbitrary obliquity, we give an expression for the total electric 
charge of the system. It does not vanish except for the perpendicular rotator. This is 
due to the often ignored point charge located at the centre of the neutron star. It is 
questionable if such solutions with huge electric charges could exist in reality except 
for configurations close to an orthogonal rotator. The charge spread over the stellar 
crust is not a tunable parameter as is often hypothesized. 
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1 INTRODUCTION 

The problem of energy losses, i.e. particle acceleration and 
radiation, from an active pulsar is closely related to the 
electric current circulation within its magnetosphere, the 
geometry of close and open magnetic field lines, the pro- 
cesses of gap formation and the way the current is closed. It 
has been argued that this current flows outside the light- 
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cylinder, where very strong dissipation arises and mani- 
fests itself as high-energy radiation. Because the electro- 
magnetic field dominates the dynamics at least close to the 
neutron star surface, low or even vanishing particle iner- 
tia and zero pressure is often assumed. This corresponds 
to a zero-th order approximate solution to the pulsar mag- 
netosphere and refered as relativistic force-free electrody- 
namics or more properly electromagnetodynamics. Steady 
axisymmetric magnetospheres have been tr eated on an ana - 
lytical basis for charge-separated plasma bv lOkamotcl (|l974l ) 
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and then exten ded to normal plasma by the same author 
|Okamotol ll975). It is well known that the force-free mag- 
netic field in the far zone well outside the light-cylinder 
should create a current sheet propagating outwards like a 
vacuum wave |Bucklevl[l977h . The force-free approxima- 
tions has been applied to magnetospheres of compact ob- 
jects, especially to those of pulsars, know n as the pul- 
sar equation, for the axisymmetric rotator |Michell Il973al : 
IScharlemann fc Wagoner|[l973;lEndeanlll974l ') with attempts 
to solve it numerically ( Gruzinov 20061 ) . An exact ana- 
lytical solution has been found for the special case of a 
magnetic monopole (|Michell Il973bl ) . Although not realis- 
tic close to the neutron star, it gives some insight into the 
far field solution, reaching asymptotically a purely radial 
structure. Many authors investigated the aligned rotator 
because axisymmetry and time independence lead to some 
drastic simplificati ons, although this does not really corre- 



spond to a pulsar ( Contopoulos et alii 19991 ; lTimokhm1l2006l ; 
^1 120061 : 1 



oimssarov 



McKinnev 20061 ). Therefore, several 



m 

thors tried to construct a force-free pulsar magnetosphere 
by direct time-dependent si mulations of Maxw ell equations 
in the general oblique case (|Spitkovskvll2006l ). Difficulties 
with outer boundary condi tions due to reflection, cured by 
perfe ctly matching layers (jKalapotharakos fc Contopoulos! 
2009) improved this algorithm by implementing some ab- 
sorbing outer boundary conditions. However, the large ratio 
between stellar to light-cylinder radius J?*/rL = 0.2 used 
in these previous works are not realistic for observed pulsar 
periods (this would correspond to a 1 ms pulsar). Moreover, 
the system has to reach a stationary state after some re- 
laxation time that needs several rotation periods and which 
is difficult to control with boundary conditions that are not 
strictly non-reflecting. However, this force- free approach w as 
recently extended by including resistivity (|Li et al.ll2012h . 

Finding analytical or semi-analytical sol utions is pos- 
sible only in a very few cases. For instance, iMestel et al.l 
l|l999l ) solved the equations for the perpendicular rotator 
without aligned currents. Therefore the linearity of the prob- 
lem allowed them to employ Fourier techniques for expand- 
ing the solution. 

However, more real istic models sho u ld lea ve the force- 
free description. Indeed fkouma fc Oogil (|2009h presented a 
two-fluid solution to the aligned rotator. Nevertheless, al- 
ternative models with non-neutral fully charged separated 
magnetospheres containing huge vacuum gaps have also 
been proposed and known as electrospheric models. The 
early simulations of iKrause- Polstorff fc Michell dl985|)_ for 
the aligned rotator have been improve d by Petri et alj (I2002T) 
and extended to oblique rotators (|McDonald fc Sheared 
l2009h . 

Back to the force-free simulations discussed above, some 
drawbacks should be pointed out. First, using a Cartesian 
coordinate system renders it painful to impose proper and 
convincing inner boundary conditions on the stellar surface, 
to a good approximation depicted as a perfect sphere. In- 
vestigating the polar cap configuration becomes therefore 
a difficult task. On the outer boundary, reflecting outgoing 
waves can hardly be avoided as explained. Therefore, re- 
laxation to a stationary state becomes difficult to recognize. 
Probably the strongest flaw is that the full expression for the 
force-free current density is not taken into account, the part 
along the magnetic field line being dropped because of its 



intricated expression including spatial derivatives, tr i cky to 
handle with finite difference schemes. IParfrev et al.l (|201ll ) 
improved the situation by implementing a pseudo-spectral 
code in spherical geometry for the axisymmetric pulsar mag- 
netosphere. Unfortunately, the code presented in this paper 
does not strictly, or at least numerically, satisfy the diver- 
genceless of the magnetic field, divB = 0, for instance by 
an explicit imposition of this condition. Spurious magnetic 
monopoles could locally be present when evolving Maxwell 
equations in time. Moreover it seems that regularity con- 
ditions at the poles were overlooked and are not ensured. 
This is probably due to inadequate expansion techniques of 
a vector field viewed as a set of scalar fields for each of its 
components. 

On a more fundamental side, lUchidal I 1997a) pre- 
sented a formalism introducing Euler potential to solve the 
time-dependent force-free equa t ions wit h speci al emphasize 
to symmetries (|Uchidalfl997bh . IPunslvl (|2003h studied the 
waves allowed by the force-free approximation and found 
that the Alfven and fast modes behave similarly to MHD 
waves in the limit of vanishing particle inertia. 

In this paper, we demonstrate how to overcome many 
of the aforementioned limitations by a proper and efficient 
expansion of a general vector field onto vector spherical 
harmonics, ensuring regularity and smoothness across the 
poles. We would like to stress that the simulations pre- 
sented here were performed on a single processor for several 
hours to days or weeks and relatively modest grids of reso- 
lution N r X JV# X N v = 257 x 32 x 64 at most. This is in 
clear contrast with finite difference simulations employing 
grids as huge as 1000 3 therefore inevitably requiring paral- 
lelization and/or supercomputers. We present full force-free 
pulsar magnetosphere solutions in nearly stationary state 
using spherical polar coordinates. The outline of this pa- 
per is as follows. The time-dependent force-free model is 
exposed and compared to a situation where the regime is 
strictly stationary, fj2] The pseudo-spectral method of solu- 
tion is presented in Sj3] The code is then tested and vali- 
dated against exact analytical solutions, SQ] Next, the mag- 
netospheric geometry and properties are described in depth 
in iJS] providing the link to the base of the striped wind, 
region believed to produce the high - energy puls ed emission 
l|Petrill2009l ; iBai fc Spitkovskvllioiol ; IPetr illioTll ) . The mag- 
netic topology is presented for an aligned, an oblique, and 
an orthogonal rotator, energy losses versus obliquity and to- 
tal electric charge. Conclusions and possible extensions are 
drawn in SJS] 



2 THE MODEL 

In this section, we briefly recall the time-dependent Maxwell 
equations in a force-free regime. Next we address the prob- 
lem of stationary force-free equations for an oblique rota- 
tor and demonstrate that the magnetic field is the main 
unknown from which all other meaningful quantities are im- 
mediately computed. Actually the problem is reduced to two 
kind of scalar fields by a vector spherical harmonic expan- 
sion. The equations for the stationary state serve as an a 
posteriori check for the time-dependent simulations to their 
closeness to stationarity. 
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2.1 Time-dependent Maxwell equations 

The time-dependent Maxwell equations in a flat space-time 
read 



d t B = 



rotE 



d t E = c 2 rotB — 



(1) 
(2) 



supplemented with the two initial conditions 

divB = (3) 



divE = 



El 
so 



(4) 



In a magnetically dominated relativistic flow such as those 
expected in pulsar magnetospheres, the force-free approxi- 
mation implies negligible particle inertia as well as negligible 
thermal pressure. The Lorentz force acting on a plasma fluid 
element must therefore vanish, i.e. 



p c E+j AB = 



(•») 



the charge density being deduced from Maxwell-Gauss equa- 
tion as 



£o div E 



(6) 



This defines the particle density number with respect to the 
electric field and not conversely. Thus we assume that there 
is a source of charge if necessary to maintain the force-free 
electromagnetic field. 

From the force- free condition and the set of Maxwel l 
equations, the current density is found to be l|Gruzinovl2006h 



E A B B 

J = PC + — 



rot B/fio — £o E ■ rot E 
IP 



B 



(7) 



The first term represents the contribution from convection of 
charges due to the electric drift, whereas the second contri- 
bution represents the current sustained along magnetic field 
lines. This expression holds only for magnetically dominated 
flows. Regions where this condition would not be fulfilled 
could in principle form behind the light-cylinder and should 
be avoid. It is difficult to design an elliptic solver passing 
through the critical points of the flow and meanwhile im- 
posing the conditions for the force-free approximation. We 
tried different iterative algorithms with several prescriptions 
behind the light- cylinder but did not get highly accurate so- 
lutions as we would expect from spectral methods. We there- 
fore decided to switch back to a time dependent formulation 
of the problem as exposed in the previous lines. These time- 
dependent Maxwell equations {J} and ((2} are integrated with 
standard explicit schemes as will be explained in section [3] 



2.2 Towards the stationary state 

It is difficult to decide whether or not the system has reached 
a steady state at the end of the run. This is especially true 
for non axisymmetric magnetospheres for which all fields 
remain time-dependent in the observer frame. It is there- 
fore useful to get a more quantitative criterion for station- 
ary as described in this paragraph. In a quasi-stationary 
state, the neutron star drags its magnetosphere and its 
wind into a rigidly corotating body with constant angu- 
lar speed Q. Mathematically speaking, this is written as 



an equivalence between partial time derivative and par- 
tial azimuthal derivative su ch that for any vector field V 
l|Muslimov fc Hardindl2005h 

9 t V = -na v V + nAV = rot(/3AV)-/3divV (8) 

where we introduced spherical polar coordinates denoted by 
(t, r, ip) and /3 = fiAr is the corotation velocity. Note that 
the partial derivative with respect to ip 



d lp {y T er + V*Z» + V' t 'e lp ) 



(9) 



includes the orthonormal basis decomposition (e T ,e^,e v ). 
Therefore this basis will also suffer some transformations 
during the differentiation process. This explains the pres- 
ence of the cross product ft A V in the second term in the 
middle part of equation (JSJ. Stated more simply, it means 
that we only need to take partial derivatives of the compo- 
nents of the vector V and not of the vector itself (both being 
different due to the curvilinear coordinate system). Thanks 
to this stationarity assumption, Maxwell- Ampere equation 
is integrated into 



-(5AB-V* 



(10) 



^ being the electric potential as measured in the stellar coro- 
tating frame. It is well known that force-free electrodynam- 
ics does not allow particle acceleration along magnetic field 
lines because E-B^B-V^^O. This also implies that 
magnetic field lines are equipotentials for ^. Moreover, re- 
calling that in the stellar interior, assumed to be a perfect 
conductor, we have 



E = —ft A B 



(11) 



implying V = 0. If there is no gap within the magneto- 
sphere, this assertion must hold in whole space. Therefore 
Eq. (|f f | remains valid everywhere, in the neutron star as 
well as in its magnetosphere. As a consequence, the elec- 
tric field is just an unessential auxiliary unknown. More 
importantly, this expression for the electric field automat- 
ically satisfies the force-free regime, i.e. E • B = 0. More- 
over this implies that at least within the light- cylinder, the 
flow remains magnetically dominated i.e. E < cB. To avoid 
unphysical solutions it is compulsory to add a significant 
toroidal component of the magnetic field outside the light- 
cylinder in order to reduce the electric drift speed. Nothing 
forbids us to extend the investigation well behind the light- 
cylinder but we were not able so far to find an efficient and 
self-consistent way to enforce this condition E < cB for the 
boundary value problem implied by the stationary magneto- 
sphere equations. On a more physical ground, in this region, 
we must either introduce an effective dissipation mechanism 
or allow for particle acceleration, abandoning the force-free 
assumption. A numerical trick to enforce magnetically dom- 
inated flows everywhere consists to reset the electric field 
to values less than cB in appropriate regions. We are then 
naturally led to a kind of relaxation process that can be ef- 
ficiently cast into a time-dependent problem. The algorithm 
would therefore be very similar to the true time-dependent 
Maxwell equations so we decided eventually to switch back 
to time-dependent simulations of the magnetosphere. 

It is readily seen that the charge density Eq.©, the 
current density, Eq.© and the electric field Eq. (|I f [) are all 
derived from the magnetic field structure. The latter is the 
only relevant unknown in a stationary force-free problem. 
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In order to investigate the closeness of the time- 
dependent solution to a stationary state, we check a pos- 
teriori that E ~ — /3AB. Moreover, the stationary Maxwell- 
Ampere equation reads 



rot B = fi j 



A E 



(12) 



which furnishes the second criterion to check the relevance 
of relaxation to a stationary state. 



3 NUMERICAL ALGORITHM 

Finding the structure of the stationary pulsar magneto- 
sphere is equivalent to relaxing the time-dependent Maxwell 
equations towards a stationary state with vanishing partial 
time derivatives d±E = <9tB = in the corotating frame or 
satisfying Eq.([8]) in the lab frame. To this end, we developed 
a pseudo-spectral collocation method in space to compute 
spatial derivatives and augmented by a standard explicit in 
time ODE integration scheme. The real strength of our code 
is the use of spherical geometry without coordinate singular- 
ity along the polar axis and proper boundary conditions on 
the stellar surface thanks to the vector spherical harmonic 
expansion exposed in many details in the appendixfX] Regu- 
larity and smoothness conditions are also automatically sat- 
isfied by the vector fields. 



3.1 Vector expansion 

A clever expansion of these vector fields is the heart of the 
code. Indeed, electric and magnetic fields are expanded into 
vector spherical harmonics (VSH) according to 

oo I 

E = Yl E (CY !m + < ) * im + ^ ) * lra ) (13) 

1 = m= — l 
oo I 



B 



J2 E ( s r m Y !ro + * !m + B« * im ) (14) 



Note that the series expansions contain a finite number of 
terms, each of them being smooth. As a consequence, the as- 
sociated vector fields will also remain smooth in the whole 
computational domain. By definition, it is impossible to rig- 
orously catch a discontinuity with spectral methods. The 
best we can do is to introduce artificially a transition layer 
of negligible thickness. From the expansion coefficients, the 
linear differential operators like div and rot can be easily 
computed in the coefficient space and then inverted back to 
real space. 

It is a special case of spectral methods, known to pos- 
sess intrinsically very low numerical dissipation and able 
to resolve sharp boundary layers (Boyd 200]]). However, if 
discontinuities arise in the solution, due to the Gibbs phe- 
nomenon, the convergence rate is drastically altered. This 
is cured by applying some filtering process to avoid aliasing 
and to smear the solution, which is equivalent to some very 
low d issipation and known as super spectral viscosity l|Mal 
1 19981 ). Vari ous kind of fi l ters c an be used as described for 
instance in lCanuto et al.l l|2006l ). 

Compared to previous works, our approach has sever- 
als advantages. First, because the expression for the current 
density is quite intricated, the field aligned component is 



usually dropped and replaced by a prescription for dissipa- 
tion. Here we implement the full expression, convection and 
field aligned contribution, reducing numerical dissipation to 
the lowest possible value. Second, outgoing wave boundary 
conditions are handled exactly by a characteristic compat- 
ibility method. Third, inner boundary conditions are pre- 
scribed on the stellar surface, i.e on a sphere. We impose 
continuous electric field components {E^,E V } as well as a 
continuous B r component. Thus the system is well posed 
and reduces to Deutsch solution for a magnetosphere with- 
out current. We now give some details about the algorithm. 



3.2 Maintaining the force- free requirements and 
the divergenceless of B 

During the time evolution of the electromagnetic field, the 
force-free conditions will be violated soon or later. It is there- 
fore necessary to reinforce the E • B = and E < cB con- 
ditions regularly. Consequently, at each time step, in a first 
stage, we subtract the magnetic field aligned electric field 
component by adjusting the new electric field E' such that 



E' = E - 



E B 



B 



(15) 



This new electric field need not satisfy the requirement 
E' < cB. In a second stage, we therefore reduce the elec- 
tric field E' by a factor E'/cB to finally get the corrected 
electric field as 



E c 



= E' 



c 2 B 2 
E' 2 



(16) 



The electric field is updated with this corrected value where 
necessary. Next, to insure a divergenceless magnetic field, 
at each time step, we project the vector B onto the sub- 
space defined by divB = simply by applying a forward 
transform to the coefficients {/,^ (r, t), g^ir, t)} according 
to Eq. ()A29[) followed by a backward transform to the vec- 
tor field B, which is analytically divergenceless by defini- 
tion. This procedure removes the longitudinal part of the 
magnetic field. 



3.3 Exact boundary conditions 

The spherical geometry of our code allows an exact enforce- 
ment of the boundary conditions on the star. Moreover, 
(pseudo)-spectral methods reflect perfectly the underlying 
analytical problem and therefore only requires the same 
boundary conditions as the original mathematical problem. 
There is no need for extra boundary conditions as would be 
the case for finite difference methods. This would make the 
problem ill-posed and the code unstable. Finite difference al- 
gorithms being much more dissipative, they will damp these 
unstable solutions! 

The correct jump conditions at the inner boundary are, 
continuity of the normal component of the magnetic field B r 
and continuity of the tangential component of the electric 
field {E$,E,p}. Explicitly, we have 



E v {t, Ri,$,<p) 



B°(t,^,<p) 

-R 1 Q sintf Br{t,d,ip) 




(17) 
(18) 
(19) 
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Note that the continuity of B r supplemented with the con- 
dition Eq. automatically implies the right boundary 
treatment of the electric field. B®(t, tp) represents the pos- 
sibly time-dependent radial magnetic field imposed by the 
star (monopole, split monopole, aligned or oblique dipole). 

As a corollary of these boundary conditions, it is impos- 
sible and even not self-consistent to force the surface charge 
density to vanish on the neutron star crust. This quantity 
will be derived a posteriori from the results of the simula- 
tions. 

The outer boundary conditions are more subtle to han- 
dle. Ideally, we would like outgoing wave conditions in order 
to prevent reflection from the artificial outer boundary. The 
appropriate technique is called Char acteristic Compatibi l- 
ity Method (CCM) and described in ICanuto et all |2007l) . 
We have implemented this technique in spherical geometry. 
Indeed, the exact radially propagating characteristics are 
known and given by 



(20) 



within an unimportant factor r. In order to forbid ingoing 
wave we ensure 



Ein 



cB v = 
■cB fl = 



whereas the other two characteristics are found by 



Ev + cB v = 
E v — cBj = 



pPDE 



cB, 



PDE 



(21) 
(22) 

(23) 
(24) 

the superscript PDE denoting the values of the electromag- 
netic field obtained by straightforward time advancing with- 
out care of any boundary condition. The new corrected val- 
ues are deduced from Eq. ((2T|) - ((24)) . We expect to reach sta- 
tionarity in the whole computational domain within a rea- 
sonable time span. 



3.4 Dealiasing and filtering 

The source term of Maxwell equations represented by the 
force-free current density Eq. introduces a strong non- 
linearity in the problem and therefore aliasing effect for 
unresolved spatial scales. In principle, the numerator be- 
ing a cubic product of unkn own functions could be cured 
by zero-padding techniques l|Canuto et al.l [20061 ). but the 
denominator renders this task impossible. Dealiasing of 
a quotient of two functions can efficiently be achieved 
by inverting a linear sy stem for one dimensional problem 
|Debusschere et aLll200 5). Unfortunately, this is intractable 
for our three-dimensional spherical problem. We prefer to 
use filtering by smoothing the electromagnetic field at each 
time step after reinforcing the force-free conditions. We use 
an eighth order exponential filter in all directions with 



a{rf) 



(25) 



r) ranges between and 1. For instance, in the radial co- 
ordinate r\ = k/(N r — 1) for k £ [O..N r — 1], k being 
the index of the coefficient Ck in the Chebyshev expansion 
~ X^fc^o" 1 Ch Tk(x)- We tried several parameter sets, a 
good compromise between accuracy and stability is a — 10 
and ,9 = 8. 



3.5 Time integration 

Pseudo-spectral methods aim at replacing a set of partial 
differential equations (PDE) by a set of ordinary differen- 
tial equations (ODE) for the unknown collocation points or 
spectral coefficients. Schematically, we write it as 



<lu 



f(t,u) 



(26) 



with appropriate initial and boundary conditions, u repre- 
sents the vector of unknown functions either evaluated at 
the collocation points or simply the spectral coefficients. 

Several methods for integrating the time evolution of 
those ODE are at hand. It is worth mentioning implicit 
versus explicit schemes and single versus multi-step march- 
ing methods. Runge-Kutta schemes remain very popular as 
a single step scheme, but their main flaw is the evalua- 
tion of the RHS very often. For spectral methods, Adam- 
Bashforth (explicit scheme) and Adam-Moulton (implicit 
scheme), known as multi-step integrators, are also widely 
employed. Their advantage is that no extra cost is needed 
to compute the right hand side. 

Although many time marching algorithms could be jus- 
tified, we decided to use a third order Adam-Bashforth 
scheme advancing a vector of unknown functions u as 



71 + 1 77 

u = u 



At ( — f n - Hf™" 1 + Af™- 2 
1 12 12 12 



(27) 



where At is the fixed time step subject to some stability 
restrictions, u n = u(nAi) and f* 1 = f (n At, u(n At)). 

Let us say a few words about the time step restric- 
tion. For the time-marching scheme to remain stable, At 
should remain in the stability region defined by the inequal- 
ity «max At ^ a Axmln where w max is the maximum speed 
of the waves, here equal to the speed of light c, a a con- 
stant close to unity depending on the integration scheme 
and Axmin the minimum grid spacing between two points. 
In our spherical computation domain, the severe constraint 
comes from the inner part because 



Ax n 



Min( Ar, r Ai9, r sin <d Aip) 



(28) 



It is usually claimed that the non-uniform grid introduced 
by the Chebyshev collocation points, which are very dense 
at the boundaries, are the most critical conditions be- 
cause Ar = 0(N^ 2 ). While this is true for Cartesian co- 
ordinates, we find actually that the grid spacing along the 
polar axis imposes an even stronger limit on the time step 
because of the presence of the factor r sin $ in front of Aip, 
depending on the resolution adopted in the radial and trans- 
verse directions. This explains why we do not tried a map- 
ping of the Chebyshe v collocation points by the popular 
Kosloff/Tal-Ezer grid (|Kosloff fc Tal-Ezer|[l993l ). Such map- 
ping unfortunately kills the spectral convergence properties 
of the meth ods, degrading it only to a second order scheme 
(|Bovdll200lh . 



4 TEST 

Before handling the general oblique rotator, we check our 
algorithm against some well known analytical solutions. The 
starting point is the Deutsch vacuum field solution for a 
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perpendicular rotator. Then, we investigate the monopole 
solutions with a central monopole magnetic field. 

In the remaining of the paper, we adopt the following 
normalization: the magnetic field strength at the light cylin- 
der in the equatorial plane is set to unity, Bl = 1, as well as 
the stellar angular velocity and the speed of light Q = c = 1, 
therefore tl = 1. 



4.1 The Deutsch solution 

The vacuum electromagnetic field of an oblique rota tor is 
known analytically since the work by iDeutschl ([1955:). Be- 
cause Maxwell equations are linear in vacuum space, solu- 
tions are found through complex analysis. Assuming a star 
of radius R* and surface magnetic field of B* with rotation 
speed f2„ and obliquity \, the complex electromagnetic field 
is easily computed and given in the appendix iBl 

We start the simulation with the perpendicular static 
dipole magnetic field (x = 90°) and zero electric field out- 
side the star, except for the crust where we enforce the inner 
boundary condition with corotating electric field, see Sec. [3] 
We performed simulations with different radial extensions of 
the computational domain and tried several resolutions. For 
the largest radial dimensions, finer grids in radius are re- 
quired, therefore increasing the number of radial collocation 
points N r . Typically, for a domain r £ [0.2, 2], a resolution 
of N r x N# x N v — 33 x 16 x 32 is sufficient to get accurate 
solutions. For the largest runs with r G [0.1, 10] we found 
that a minimum resolution of N r x N# x N v = 257 x 16 x 32 
was necessary. 

We let the system evolve for two rotational period of 
the pulsar. The final magnetic field line configuration in the 
equatorial plane is shown as dashed blues curves in fig. [T] 
and compared to the analytical solution depicted as solid 
red curves. The numerical solution is very close to the exact 
solution everywhere, they are hardly distinguishable. The 
associated Poynting flux is shown in fig. [5] It is constant 
and equal to the Deutsch value everywhere except close to 
the star. Nevertheless after a thin boundary layer, the flux 
relaxes sharply to the exact stationary value of Deutsch. 
In the outer parts of the domain, the system evolved to a 
nearly steady state configuration and no significant reflec- 
tions have to be reported. The Poynting flux is normalized 
to the monopole solution flux, see Eq. (|32|l . It is a bit less 
than unity because it depends on the ratio i?*/rL and tends 
to Eq. (|32p when R t /r^ tends to zero. 

This first test demonstrates the ability of our code to 
handle non axisymmetric configuration with high accuracy 
and to relax the system to an almost stationary state in a 
vacuum region. 



Deutsch field lines for the orthogonal rotator 




Figure 1. Magnetic field lines of the perpendicular Deutsch field 
in the equatorial plane. The exact analytical solution (solid red 
lines) is compared to the time-dependent simulation (dashed blue 
lines). The overlapping is excellent in the whole region where the 
relaxation to a quasi-stationary state is nearly achieved. 
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Figure 2. Poynting flux across the sphere of radius r. The com- 
puted flux (blue line) is almost constant and equal to its analyt- 
ical value (red line). The solution settled down to the stationary 
Deutsch field even close to the outer boundary. 



In terms of a vector spherical harmonic (VSH) expansion, 
this magnetic field is expressed as 



B 



r 2 



•0io I 



(30) 



where 



4.2 The monopole solution 

Next we tackle the problem of an axisymmetric fo rce- free 
flow kn own as the monopole field introduced by iMichell 
(1973b). We recall that this monopole solution is given by 



B 



D r L . q -* 

i3L — sm w e ¥ 



(29) 



ffioM = \HtBl — 
V 3 r 



(31) 



all other coefficients being equal to zero. See the appendix |A"1 
for a detailed exposure of the VSH, their definition and use- 
ful properties. 

The run is started with a pure monopolar magnetic 
field and zero electric field outside the star, as before. 
An Alfven wave is launched from the stellar surface and 
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Figure 3. The relevant VSH coefficient gf (r) for the monopole 
magnetic field expansion in the domain r/r-^ S [0.1, 10]. For ease 
to compare, we applied a normalization factor of rj ^/8tt/3. All 
other coefficients gfLAr) with the same normalization factor are 
less than 10~ 6 . Results are shown for a resolution of N r X N$ X 



257 x 16 x 32. 
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Figure 4. Poynting flux across the simulation sphere correspond- 
ing to the numerical solution of fig. [3j It is approximately constant 
as expected from energy conservation arguments. The same res- 
olution as in fig. \3\ applies. 



quickly forces the system to a stationary state. Here also, 
as expected, larger radial extension implies finer grid in ra- 
dius. Typically, for a domain r G [0.2,2], a resolution of 
N r x N# x N v — 33 x 16 x 32 is sufficient to get accurate so- 
lutions. For the largest runs with r £ [0.1, 10] we found that 
a minimum resolution of N r x N& x N v = 257 x 16 x 32 was 
necessary. The relevant coefficient gf (r) is shown in fig. [3] 
for comparison with the analytical solution. The numerical 
solution is very close to the exact solution in the whole com- 
putational box. The associated Poynting flux, normalized to 
the exact monopole value of 

■ fi* B L r L (32) 



■ fit bI r\ = 



3 fio c 3 * * * 3 fio c 3 

is depicted in fig. [4] It is almost constant and equal to the 
monopole value as expected. Note that there is no spuri- 
ous reflection in the outer regions, we are very close to a 
steady state configuration. We tried higher resolutions for 
instance N r x N$ x N v = 257 x 32 x 64 without significant 
improvement of the computed solutions. All spatial scales 
have been resolved with the coarsest grid. That such a low 
resolution is sufficient is no surprise because the monopole 



solution does only contain two coefficients: the monopole ra- 
dial part foo(r) = -Bl r^/r 2 and the toroidal part gu>{r), so 
there exist at most one mode 1 = 1. 



5 RESULTS 

We next look at the oblique pulsar magnetosphere, espe- 
cially the aligned rotator, studied several times in the past 
by several authors. New results about the perpendicular ro- 
tator and its link to the striped wind will also be shown. 



5.1 Simulation setup 

The rotating neutron star has a radius R\. The computa- 
tional domain extends from this inner boundary R± to an 
outer boundary shell of radius R2 > Ri- Note that the star 
does not belong to the simulation box except for its surface 
in order to fix the inner boundary conditions. For concrete- 
ness, we start from Ri = 0.2 rL corresponding to a 1 ms 
pulsar down to R = 0.1 tl describing a 2 ms pulsar. Like- 
wise the extension of the simulation sphere is variable and 
equal to R2 = 2 tl up to R2 = 10 tl to investigate the 
beginning of the striped wind. The inclination angle of the 
magnetic moment with respect to the rotation axis was set 
to x= {0°, 30°, 60°, 90°}. 

We performed several simulation sets with different fil- 
ters, the most common being Cesaro, r aised cosinus and 2nd , 
4th, 8th-order exponential filter (see ICanuto et all |2006)) 
as well as different spatial resolutions. We found that for 
spatially resolved discretization, the algorithm converges 
quickly to a stationary state. We emphasize that the so- 
lutions become independent of the choice of a filter as it 
should be, given us confidence about the convergence and 
reliability of our results. 



5.2 Magnetospheric structure 

First, let us have a look on the magnetic field geometry. Let 
us start with the special case of an aligned rotator. The 3D 
magnetic field lines are shown in fig. [5] The inner boundary 
is at Ri = 0.2 tl and the outer boundary at R2 = 2j"l. In- 
side the light-cylinder, the geometry is close to the static 
dipole without significant toroidal magnetic field compo- 
nent. Along the polar caps and outside the light-cylinder, 
we recognize a monopole structure with a significant toroidal 
component. The aligned rotator, although being simple be- 
cause axisymmetric and therefore more tractable analyt- 
ically and computationally, is of no relevance for pulsar 
theory because it is not able to reproduce the pulsar phe- 
nomenon. Inclined models are true pulsars. Thus, we give an 
example of oblique rotator as shown in fig. [6] with x — 60°. 
The weakly disturbed static dipole regime inside the light- 
cylinder contrasts with the open field lines passing through 
the light- cylinder as is costumary for the force-free magneto- 
sphere. Visualization and exploiting data is cumbersome for 
the full 3D geometry, so we do not discuss further this pecu- 
liar case. The orthogonal rotator with its inherent symmetry 
about the equatorial plane represents the best configuration 
to study a realistic pulsar magnetosphere. We show such ro- 
tator for which vizualisation is easier because some field lines 
are entirely contained in the equatorial plane. These lines 
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Figure 5. 3D magnetic field lines for the aligned rotator. The 
inner and outer spherical boundaries are located respectively at 
R\ = 0.2 r L and R 2 = 2r L . 



z/rL 




Figure 6. 3D magnetic field lines for the inclined rotator with 
X = 60°. The inner and outer spherical boundaries are located 
respectively at Ri = 0.2 rL and Ri = 2rjj. 

are shown in fig. [7] The inner boundary is at Ri — 0.2 tl 
and the outer boundary at Ri — 5tl- The magnetic topol- 
ogy is reminiscent of the Deutsch field solution with some 
regions of strong deviation due to the magnetospheric cur- 
rent. Inside the light-cylinder (black circle), the geometry 
is close to the static dipole whereas outside, the launch of 
the striped wind is seen with its spiral structure in regions 
of strong magnetic field gradient. In order to compare the 
orthogonal force-free field with the vacuum field, we did a 
large scale simulation with Ri = 0.1 rL and Ri = 10 tl and 
shown in fig. [8] The shape of the spiral current sheet in the 
force-free striped wind follows a geometric shape reminis- 
cent to the Deutsch solution. The locus of the current sheet 
in the striped wind according to the rotational phase of the 
neutron star has strong implications on the association be- 



Equatorial magnetic field lines for the orthogonal rotator 
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x 

Figure 7. Equatorial magnetic field lines for the perpendicular 
rotator. The light cylinder is shown as a black circle. The inner 
and outer spherical boundaries are located respectively at R\ = 
0.2 rL and R2 = 5r^. 



Equatorial magnetic field lines for the orthogonal rotator 
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Figure 8. Equatorial magnetic field lines for the perpendicu- 
lar force- free rotator (red solid line), compared to the vacuum 
Deutsch field (dashed blue line). The light cylinder is shown as a 
black circle. The inner and outer spherical boundaries are located 
respectively at Ri = 0.1 rL and Ri = 10 rL- Note the similarity 
between both structures. 



tween radio and high energy emission. The time lag between 
the radio pulses and the gamma-ray pulses a re direct ob - 
servational consequences and can be checked (|Petrill201ll ). 
Inspecting the Poynting flux dependence with respect to ra- 
dius, fig. |9j we observed a tendency to dissipate the electro- 
magnetic flux well outside the light-cylinder. The employed 
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Figure 9. The Poynting flux across the whole simulation sphere 
for the orthogonal rotator of fig. \E\ The Poynting flux is nearly 
constant but artificial viscosity due to filtering slightly alters en- 
ergy conservation outside the light-cylinder. 
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Figure 10. The Poynting flux across the whole simulation sphere 
as a function of the obliquity \. The flux is nearly constant for 
the perpendicular rotator. Artificial dissipation due to filtering 
becomes significant for almost aligned rotators outside the light- 
cylinder. 



resolution should still be increased but must wait more com- 
putational power. 

5.3 Energy losses and electric charge 

Another important issue is the spin-down energy of the pul- 
sar by magnetic braking. Thanks to the knowledge of the 
magnetospheric structure, we can investigate the energet- 
ics of the pulsar quantitatively. Current wisdom assumes 
that the Poynting flux escaping from the force-free mag- 
netosphere is of the same order of magnitude as the vac- 
uum dipole rotator. We check this assertion by plotting the 
Poynting flux across the simulation sphere as a function 
of the obliquity with a simulation sphere extending from 
Ri — 0.2rL to 7?2 = 27"l, fig. 1101 normalized with respect 
to the magneto-dipole losses for an orthogonal rotator ex- 
pressed as 

Ldip - — (33) 

Within a factor of several unity, the luminosity is the same 
for vacuum dipole and force-free magnetosphere, fig. 1101 The 
flux across the light-cylinder, i.e. before artificial viscosity 
due to filtering becomes significant can approximately be 
fitted with, fig. QT] 

L sp » 1.5L dip (l + sin 2 x) (34) 

Thus the orthogonal force-free magnetosphere radiates 
3 times more than the vacuum analog whereas the aligned 
magnetosphere radiates 50% of the vacuum perpendicular 
ro tator. These result s agree with the functional form found 
bv lSpitkovskvl (|2006T ) although we took into account the full 
expression for the current density. Surprisingly, we found 
it easier to look for the perpendicular rotator than for the 
aligned one. Indeed, in the latter case, the current sheet in 
the equatorial plane is coarsely resolved by our current grid 
resolution. According to the lower blue curve in fig. 1101 non 
negligible artificial dissipation is present and smears out the 
discontinuity. On the contrary, the strong displacement cur- 
rent in the former case leads to more gentle current sheet 
much better described by the coarse grid. Scrutiny of the 
upper green curve in fig. [10] demonstrates indeed that the 
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Figure 11. The Poynting flux across the sphere of radius n, as 
a function of the obliquity \- Simulation results are depicted by 
red triangles and the fit is shown in blue solid line. 



Poynting flux remains nearly constant with radius. Energy is 
conserved to good ac curacy as is expect ed from the force-free 
evolution equations |Komissarov|[20 0^j) . In order to circum- 
vent high dissipation rate, the only efficient remedy would 
be to increase the grid resolution. Unfortunately, for the 
moment we are unable to run such high-resolution simu- 
lation on our computer. We plan to improve our code by 
performing some optimizations and by using parallelization 
techniques. 

Finally, the total electric charge enclosed by a spher- 
ical shell centered at the pulsar location is computed by 
the electric flux across it. It is too rarely stated that the 
total charge of the neutron star and its magnetosphere is 
not equal to zero but very close to the value of the central 
point charge Q c , located in the middle of the star (see ap- 
pendix iBl for the definition of this charge). Indeed, the total 
charge of the system with respect to obliquity \ is shown 
in fig. 1121 Only in the special case of an orthogonal rotator 
is this charge exactly zero. This is true for all ratio R^/r^. 
The fit is simply given by Qg c ~ 0.9 Q c , fig- 1131 Here also, 
the dependence on cosx is reminiscent of the central point 
charge Q c . 
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MHD equations. This would bring us closer to reality in the 
striped structure. Moreover, the magnetic field originates 
from a compact object of size only twice its Schwarzschild 
radius. Therefore, we expect not only magnetospheric cur- 
rents distortion of the magnetic topology, but also general 
relativistic effects, especially important in the vicinity of the 
stellar surface and very useful for elucidating the precise ge- 
ometry of the polar caps. These perturbations due to curved 
space and dragging of inertial frames, should also be taken 
into account for an accurate investigation of the emission 
region close to the magnetic poles. 



Figure 12. Total electric charge Qtot in the simulation sphere as 
a function of the obliquity \ an d normalized against the central 
point charge for the aligned rotator Q c / cos \. Simulation param- 
eters are exactly the same as in fig. 1101 
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Figure 13. Total electric charge across the sphere of radius tl 
as a function of the obliquity \- Simulation results are depicted 
by red triangles and the fit is shown in blue solid line. 

6 CONCLUSION 

We solved the full set of time-dependent force-free equations 
for a rotating oblique pulsar magnetosphere. We tested our 
pseudo-spectral collocation code and demonstrated that ac- 
curate numerical solutions can be found for smooth fields. 
Interestingly, we found it easier to compute the solution for 
the perpendicular rotator because the flaw of any spectral 
method to catch discontinuities is harmless in that particular 
geometry. It seems that the presence of a strong displace- 
ment current facilitates the computation in the vicinity of 
the current sheet. Unfortunately, the current version of our 
code does not allow to perform very high resolution on a 
single processor in a reasonable time to reduce the Gibbs 
phenomenon for non-smooth configurations as those close 
to the aligned rotator. We plan to improve the code by op- 
timization and parallelization techniques in the near future. 
The boundary conditions could certainly also benefit from 
some more elaborated methods like a sponge layer and non 
alteration of the boundary value by the filtering process (not 
yet implemented). 

We think that the most physical satisfactory way to 
get rid of the Gibbs phenomenon related to the thin current 
sheet would be to alleviate the ideal MHD approximation by 
including particle inertia. Our code is sufficiently versatile 
and can be straightforwardly extended to solve the full set of 
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APPENDIX A: SCALAR AND VECTOR SPHERICAL HARMONICS 

In this appendix, we expose in some details our pseudo-spectral algorithm to compute approximate analytical solutions to 
the force-free pulsar magnetosphere according to the time-dependent Maxwell equations. 

The heart of our scheme resides in the vector spherical harmonic expansion of the two unknown vector fields, the electric 
and magnetic fields. In the following paragraphs, we discuss in depth the properties of this vector basis. But let us first quickly 
remind the scalar spherical harmonics. 



Al Scalar spherical harmonics 

The scalar spherical harmonics are defined according to the associated Legendre functions P/ 71 by 



Y lm{ ») = ^ {^^(cos,) (Af) 

y, m (0,p) = Y im (tf)e lm¥> (A2) 

The square root factor imposes normalization of these functions. They are very valuable functions because eigenfunctions of 
the Laplacian operator in spherical coordinates such that 

AYi m (#, v) = - 1 ( * r t X) P»m( J, y) (A3) 

A useful symmetry property is 

if) = (-i) m yf m (tf, V ) (A4) 

Note that the Yi m form an orthonormal basis on the two-dimensional sphere. Any smooth and regular scalar field can be 
expanded on this sphere thanks to the spherical harmonics. As explained in the next paragraph, this statement can be 
generalized to a vector field defined on the sphere. 

For completeness, we give the first few Yj m up to m = 3 

Y 00 ' 



20F 



Yl0 = yicosw 



Y 20 = -\ - (-l + 3cos[tf] 2 ) 
4 y 7r 

Y 21 = -i e ^ A /|^co S [i?] S in[i?] 

Z V 27T 
v 1 2i<p / 15 . r Q1 2 

22 = 4 6 V2^ Sm[ 1 
Y 30 = i^(-3cos[^] + 5cos[^] 3 ) 

Y31 = -i e ^y^(-l + 5cos[i9] 2 )sin[' 

Y 32 = ie 2l ^/^cos[tf]sin[#] 2 

4 V 27T 

1 3i» /35 3 



Y 33 = --e^/^rint* 



A2 Vector spherical harmonics 

Vector spherical harmonics (VSH) are often used in atomic physics to compute electronic level in atoms. Quantum physicists 
therefore prefer to work with a normalization of the vectors involving the complex number i. Here we adopt a slightly different 
definition of these vectors by removing this unessential constant factor i. More precisely, the three sets of vector spherical 
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harmonics we use are defined by 



r 



A VYj m 



(A5) 
(A6) 

(A7) 



a +i) 

Any smooth three-dimensional vector field E admits an expansion onto these vectors according to 

oo I 

E(r,0,<p) =£ X) (Wta + S«(r)*, m + ^(r)# im ) (A8) 

(=0 m=-( 

In order to avoid complex numbers, it is sometimes useful to employ instead the trigonometric representation like for instance 

r c/s 



,m V sin mtpj 

and similar expressions for , vf^y. 

To link our definition to other authors, note that in quantum mechanics the <I>; m are named X; m and defined by 



(A9) 



X, 



v/Tp+T) 



r A VF; m = -i * iri 



(A10) 



A3 Properties 

The vector spherical harmonics share some useful properties with respect to their spatial derivatives. First, assume a 3D scalar 
field expanded onto the scalar spherical harmonics such that 



1=0 m=—l 

Then, its gradient expanded onto the VSH becomes 



oc I / 

w = XE 

1=0 m=—l \ 



— - I|m 1 <plm Vim 

or r 



The action of the same gradient on the VSH gives the divergence of any vector field E as 

oo I 



1=0 m = -l 



{r 2 E r lm ) 



and for the curl 

oc I 

VAE = EE 

1=0 m=-i 

For each component taken individually, we find for the divergence 



V-(/(r)Y im ) = 



V • (/(r) 

V ■ (/W*Im) 

and for the curl 
Vx(/(r)Y lm ) 

V x (/(r)*i m ) 



v/Kl + i) 



= 



1 d 



V x (/(r)* 



/m J 



1 9 



r r dr 

Finally, define the radial differential operator T>i by 



Vi = — — (r 2 — 

r 2 dr V dr 



lm ] £m 



Z(7 + l) 



(All) 



(A12) 



(A13) 



(A14) 



(A15) 

(A16) 
(A17) 



(A18) 
(A19) 

(A20) 
(A21) 
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The VSH noted <I>; m are eigenvectors in the linear algebra meaning, of the vector Laplacian operator A. Indeed, it is straight- 
forward to show that 

A[/(r)* !m ] =£>;[/] * im (A22) 
It is thus an extension of the properties of the scalar spherical harmonics to the realm of 3D vectors. 



A4 Expansion of a vector field onto VSH 

From the above discussion, the components of any vector field can be computed according to the three underlying equalities 
E r = ^E[ m Y lm (A23) 

l,m 

div^E = £-Y5EI) £ Wy !ra (A24) 



rotter = ^-^SM E ^Y lm (A25) 

l,m 

where div^i? means taking only the angular part of the divergence. More explicitly, from the definition of the differential 
operators, we get 

E r = Y. E ^Yim (A26) 

l,m 

-L-doismdE^ + ^—apEv = Y'-y/l(l + l)E™Yi m (A27) 

l.m 

-L-d#(sm&E v )- J—d^Eo = Y-^l(l + l)E^Y lm (A28) 
smi9 smi9 ^ 

l,m 

Finding the components of E is therefore equivalent to finding the expansion coefficients of three scalar fields onto the scalar 
spherical harmonics. This procedure works for any vector field. However, the magnetic field being divergenceless, only two 
of the three components are independent. It is therefore judicious to deal properly with those kind of fields by analytically 
enforcing the condition on the divergence as explained below. 



A5 Expansion of a divergenceless vector field 

Any divergenceless vector field is efficiently developed onto the vector spherical harmonic orthonormal basis. This will be the 
case for the magnetic field in our algorithm. 

Assume that the vector field V is divergenceless. It is helpful to introduce two scalar functions fi m (r, t) and gim{r, i) such 
that the decomposition immediately implies the property of divergenceless field. This is achieved by writing 

oo I 

V(r,#,<p,t)=J2 E (rot[/ im (r,t)$ !ro ]+ ff!m (r,t)# !m ) (A29) 

£=1 m= — l 

This expression automatically and analytically enforces the condition div V = 0. Let us quickly draw the way to compute 
these functions. The transformation from the spherical components to the functions (fi m ,gim) is given by 



V-e r = f; J2 ~^±Afi m Y lm (A30) 



; = 1 m=—l 

oo I 



(rot V). el = E E -^^-girnYirn (A31) 

1 = 1 m — — l 

Thus, it is sufficient to expand again the radial component of the vector and its curl onto scalar spherical harmonics. We get 

rV r = X)-V'(i + l)/»m*Jm (A32) 

J—dtismdVv)- -J-;d V V0 = Y'-y/l(l + l)gi m Y lm (A33) 
sin w sin 17 ' 
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The functions {fim,gim} are related to the general expansion Eq. (|A8|) by 

= _V[([±i) /!m (A34) 

r 

= -;3 r (r/ im ) (A35) 

^ = Sh» (A36) 



A6 The first few VSH 

Here also, for completeness, we give the first few VSH which can be useful for analytical calculations. 
The VSH S&im are given in the orthonormal basis {e r , e#, e v } by 
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whereas the <&( m are expressed as 
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APPENDIX B: DEUTSCH FIELD 



For comparison with numerical simulations, we recall the exact complex expressions for the Deutsch vacuum field solution in 
the general oblique case. The magnetic field, dipolar close to the neutron star, is given by 
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The induced electric field is then 
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k = 1/fL is the wavenumber and /i, are the spherical Hankel functions of order I satisfying the outgoing wave conditions, see 
for instance lArfken fc Weberl (| 19951 ). The physical solution is found by taking the real parts of each component, it encompasses 
a linear combination of the vacuum aligned dipole field and the vacuum orthogonal rotator with respective weights cos \ and 
sinx- To complete the solution, we add a monopolar electric field contribution due to the stellar surface charge such that 
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where Q, is the total electric charge of the star and 
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is the central point charge inside the star. The above expressions reduce to the static oblique dipole for small distances r < tl 
(static zone). 
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